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We present a lattice study of up to = 20 unitary fermions confined to a harmonic trap. Our 
preliminary results show better than 1 % agreement with high precision solutions to the many-body 
Schrodinger equation for up to = 6. We are able to make predictions for larger which were 
inaccessible by the Hamiltonian approach due to computational limitations. Harmonic traps are 
used experimentally to study cold atoms tuned to a Feshbach resonance. We show that they also 
provide certain benefits to numerical studies of many-body correlators on the lattice. In particular, 
we anticipate that the methods described here could be used for studying nuclear physics. 
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1. Introduction 

The experimental study of trapped, ultracold atoms has provided impetus to explore such 
systems computationally, particularly near the limit of unitarity. In this regime, all scales relevant 
to the interaction vanish, corresponding to a zero range two-particle interaction with an infinite 
s-wave scattering length. This insensitivity to the scale of the interaction allows one to extract 
physical quantities which are universal in nature. Thus, in addition to being directly relevant to cold 
atom experiments, numerical studies of this system serve as a stepping stone to more complicated 
systems such as nuclei, where s-wave scattering lengths are unnaturally large compared to the range 
of the interaction [1,2]. 

A recently developed lattice method for describing unitary fermions [3] allows us to explore 
several issues of wide interest in lattice field theory: in particular, how to extract the properties 
of conformal systems from calculations at finite lattice spacing, volume, and particle density, and 
how to construct optimal interpolating fields for strongly interacting, many-fermion systems. In 
addition, we demonstrate that it is possible to obtain a very clean signal for up to = 20 fermions 
using modest computational resources. 

Our new method is presented in detail in companion proceedings [3]. In these proceedings, we 
report a preliminary study of up to = 20 unitary fermions trapped in a harmonic potential. We 
can benchmark our method and systematic errors for up to = 6 against high precision solutions 
to the many-body Schrodinger equation, achieving agreement at the 1% level. We believe this is 
the first microscopic study to explore higher without invoking a variational principle or requiring 
costly importance sampling. We will then present preliminary results for the Bertsch parameter 
(^) and the pairing gap. In another companion proceedings, we will present results for untrapped 
unitary fermions [4]. 

2. Lattice Theory 

Details of our lattice theory, along with a discussion of numerical costs and benefits, can 
be found in [3]. Our approach is a highly improved version of the theory presented in [5] at zero 
chemical potential and zero temperature. As outlined in [3], we are able to employ a simple iterative 
procedure for propagator production, whose form is 



where D is the kinetic part of the transfer matrix and X represents the four-fermion interaction, 
which is produced stochastically through an auxiliary Z2 field. This interaction is highly tuned to 
produce cot 5o ^ 1 , where 5o is the s-wave scattering phase shift, for two fermions in a finite 
box up to arbitrarily high momenta. 

We include an external potential through the addition of a field-independent term to our inter- 
action operator: 



K-\t,Qi)=D-^X{'1-1)D-^ ■■■D-^X (0)D 



(2.1) 



(2.2) 
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This introduces G(hx) errors, where hx is the temporal lattice spacing ^ We will discuss systematic 
errors in the next section. 

To build our sources, we begin with the set of single particle wavef unctions, for free 
fermions in a harmonic trap, where = {(0,0,0), (1,0,0), (0, 1,0), •• •} in the Cartesian basis. 
While a simple Slater determinant of such states is sufficient to find the ground state for the trapped 
system, a more sophisticated approach involving two-body correlations gives superior overlap. To 
build pairing into the system we use the solution for two particles at unitary in a trap, 

^-(x2+/)/(2L2) 

{x,y\\^PAiR] °^ . . , (2.3) 

V-y\ 

where Lo is the characteristic trap size, defined in the next section. In order to avoid an G{y ) 
increase in computation time, we use single particles states at the source and two particle states at 
the sink to form an {N /2) x {N/2) Slater determinant as in [6]: 

CN{T) = {det{Sij)), Sij = {YPAiR\K-\t-0)\ninj), (2.4) 

where represents the direct product \nj). For odd N, the two-body state in the highest 
energy level is replaced by the corresponding single particle state. 

3. Systematic Errors 

The relevant scales of the trapped system inside a periodic box are the box size, L, the charac- 
teristic size of the trap, Lq = (kM)^^/^, and the oscillator frequency, CO = \J k/M, where K and M 
are the spring constant and fermion mass respectively. To reach the continuum and infinite volume 
limits we require bs <^ Lq <^ L, where is the spatial lattice spacing, and b^ <C (O. Thus, our sys- 
tematic errors will be controlled by the dimensionless quantities (Ob^, bs/La, and Lq/L. To balance 
the need for small temporal discretization errors with the computational cost associated with the 
number of time steps required to reach the ground state, we have chosen (obx = 0.0136 for this 
preliminary study. 

For a given box size, the choice of Lq must take into account both discretization errors and 
finite volume errors. The expectation is that for small L^/bs spatial discretization errors will dom- 
inate. Since the discretization is implemented as a hard cutoff in momentum space, the sensitivity 
of the state to the infinite potential at the boundary results in an increase in the associated energy. 
Conversely, for large Lq / bs finite volume errors will dominate. The periodic boundary conditions 
in space result in attractive interactions from image particles, causing a decrease in energy. Thus, 
a scan of the ground state energies as a function of Li^/bs is necessary to determine at which value 
we can minimize both types of error. Fig. 1 (left) presents our findings. 

To ensure small systematic en^ors for larger N, we have chosen a value of L^/bs = 4.0, and 
performed the calculation at L/bs = 24,32,48,64 (Fig. 1, right). Based on our results for different 
volumes, we expect that for all A'^ used in this study, at L/bs = 64 the finite volume errors should 
be less than ~ 0.5%. 

'in a future publication we will adopt the formX(T) i Ei=i(''^^/2) x(T)e^ i Ei=i(''^^/-) , which has discretiza- 
tion errors of 6'(b^). 
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Figure 1 : Ground state energies as a function of Lq / bs 
at fixed L/bs = 32 (left) and L/b^ at fixed Lo/b, = 4 
(right) for various A'^. Solid lines are results from [7]. 
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4. Results 

For to N = 20 fermions, we have generated approximately 10 million configurations for 
each of the four volumes, using a total of 350 thousand CPU hours. We have found excellent 
signal for all N, with little signal-to-noise problem (Fig. 2). As discussed in [4], we find that 
results for unitary fermions in the absence of a trap are far more difficult to obtain. We believe 
that the spatial confinement enforced by the harmonic trap provides numerical benefits to lattice 
calculations of interacting systems, in particular early plateaus and far less sensitivity to the choice 
of interpolating operator. 

To extract the ground state energies, we resampled the data using the bootstrapping tech- 
nique. We then performed correlated fits to the plateau regions. Our fitting systematic errors 
were found by varying the endpoints of the fitting region by At = 2. Our results are presented 
in Fig. 3. We also show the results from two fixed-node calculations: a Green's function Monte 
Carlo (GFMC) approach [8] and a Diffusion Monte Carlo (FN-DMC) approach [9]. By using the 
fixed-node constraint along with a variational principle, both of these methods can only provide 
upper bounds on the ground state energies. We find that our energies are consistently lower than 
those obtained using both of these methods, and that this discrepancy grows with N, with N = 20 
being approximately 7% lower than the GFMC approach, and 3% lower than FN-DMC. Our values 
and their coiTcsponding statistical and systematic errors are reported in Table 1 . 

By performing correlated fits to the data for different A'^ we can extract two additional quantities 
of interest to very high precision. The first is the quantity known as the Bertsch parameter. For both 
unitary and non-interacting fermions in a box, the only relevant scale is the density, thus the ratio of 
their ground state energies must be a constant, ^ . With a trap present one can show, using the Local 
Density Approximation (LDA), that the energy of the interacting system in the thermodynamic 
hmit is given by [10]: 

Em (N, CO) = y/^Efree {N, co) . (4. 1) 
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Figure 2: Effective mass plots and fits (solid lines) for various N. The red point is the exact ground 
state energy for A'^ = 3. 



Our results as a function of are shown in Fig. 4 for L/bs = 64. One sees clear shell structure 
for the full range of A'^ considered. This indicates that we are probably not sufficiently close to 
the thermodynamic limit to extract an accurate Bertsch parameter. The inset shows con^elated 
fits to the second shell as a function of L/b^, with our preliminary value from the largest volume 
being t, = 0.450(1). For comparison, the data from GFMC and FN-DMC give t, = 0.50 and 
= 0.46, respectively. While the Bertsch parameter we extract is lower than that from variational 
parameter calculations, the value is still higher than calculations with no external potential and 
extractions from experiment. For example, a recent GFMC calculation finds = 0.40(1) [11], 
while an experimental study finds t, = 0.39(2) and = 0.41(2) using two different extraction 
methods [12]. Similarly, using our lattice method we find a preliminary value of 0.412(4) for the 
untrapped system [4]. This could indicate a slower convergence to the thermodynamic limit for 
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Table 1: Ground state energies as a function of A'^ in units of (O (preliminary). The first eiTor is 
statistical while the second is a systematic eiTor from the choice of fitting region. For N <6 v/e 
also present high-precision results from [7]. 
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Figure 3: Ground state energies as a function of N in units of CO (preliminary). For comparison, we 
show results from GFMC [8] and FN-DMC [9]. 



the trapped case, where the LDA must be invoked. In future studies, we plan to go to higher N to 
obtain a more accurate prediction for the Bertsch parameter. 
The second quantity of interest is the pairing gap, 

L{N)=E{N)-]^[E{N-\)-E{N + \)] (for N odd). (4.2) 

The gap quantifies half of the energy required to break a fermion pair. We present our results as 
a function of A'^ in Fig. 5. For comparison, we also show the results from [9]. By performing a 
con^elated fit in A'^ we are able to achieve statistical errors which are an order of magnitude smaller 
than the previous calculation. The gap shows a clear shell structure for the range of considered. 
Due to this shell structure, we will need to extend our calculation to much larger N before we are 
able verify the N'^I'^ behavior proposed in [13]. 



5. Conclusions 

We have performed a study of unitary fermions in a trap which match high precision results 
to within 1%. We have shown that the signal-to-noise problem for our lattice method is well under 
control, and predict that results for at least A'^ = 40 will be numerically feasible. In addition, we 
have performed a preUminary study of systematic errors. 

In a future publication, in addition to extending the study to larger A^, we will improve the 
discretization of our potential, introduce a Galilean invariant interaction to ensure that p cot 5o = 
for boosted pairs of particles, and perform an even more extensive analysis of systematic errors. 
Our ability to precisely tune p cot 5q suggests future applications to nuclear systems. The use of 
harmonic traps for extracting nuclear physics has also gained interest recently [14]. Due to the 
numerical benefits of including a harmonic trap, our formulation lends itself to such studies. 
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Figure A: B, = [Eint{N,(o)/Efree{N,(o)]^ as a Figure 5: The pairing gap as a function of A/^. Er- 

function of N. Error bars are smaller than data ror bars are smaller than data points. For com- 

points. Inset: Constant fit to the second shell as parison, the results from FN-DMC [9]. 
a function of L/b,,.. 



Acknowledgments 

This work was supported by U. S. Department of Energy grants DE-FG02-92ER40699 (to M. 
G. E.) and DE-FG02-00ER41 132 (to D. B. K., J-W. L. and A. N. N.). This research utiUzed re- 
sources at the New York Center for Computational Sciences at Stony Brook University/Brookhaven 
National Laboratory which is supported by the U.S. Department of Energy under Contract No. DE- 
AC02-98CH10886 and by the State of New York. The authors would like to thank D. Blume, et al. 
for granting permission to present their results before publication. 

References 

[1] D. B. Kaplan, M. J. Savage, and M. B. Wise, Nucl. Phys. B534, 329 (1998), nucl-th/9802075. 

[2] D. B. Kaplan, M. J. Savage, and M. B. Wise, Phys. Lett. B424, 390 (1998), nucl-th/9801034. 

[3] M. G. Endres, D. B. Kaplan, J.-W. Lee, and A. N. Nicholson, PoS LATTICE2010, 182 (2010). 

[4] J.-W. Lee, M. G. Endres, D. B. Kaplan, and A. N. Nicholson, PoS LATTICE2010, 197 (2010). 

[5] J.-W. Chen and D. B. Kaplan, Phys. Rev. Lett. 92, 257002 (2004), hep-lat/0308016. 

[6] J. Carlson, S. Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003). 

[7] D. Blume and K. M. Daily, ArXiv e-prints (2010), 1008.3191. 

[8] S. Y. Chang and G. R Bertsch, Phys. Rev. A76, 021603 (2007), physics/0703190. 

[9] D. Blume, J. von Stecher, and C. H. Greene, Physical Review Letters 99, 233201 (2007), 0708.2734. 
[10] T. Papenbrock, Phys. Rev. A72, 041603 (2005), ai-Xiv:cond-mat/0507183. 
[11] A. Gezerlis and J. Carlson, Phys. Rev. C77, 032801 (2008), 0711.3006. 
[12] L. Luo and J. E. Thomas, Journal of Low Temperature Physics 154, 1 (2009), 08 1 1 . 1 159. 
[13] D. T Son, ArXiv e-prints (2007), 0707.1851. 

[14] T. Luu, M. J. Savage, A. Schwenk, and J. R Vai-y, Phys. Rev. C82, 034003 (2010), 1006.0427. 



7 



